#细胞比例计算及可视化
#https://zhuanlan.zhihu.com/p/478520665

#细胞比例计算及可视化

#细胞比例####
setwd("C:/Users/forbing36/Desktop/单细胞生信/GSE212966")
scedata <- readRDS("./data/temp/Only_T_cluster_id_test.rds")
# 修改ident的顺序，解决后期作图排序
levels(x = scedata)
levels(x = scedata) <- c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_CTL", "CD8_Trm", "CD8_Te")
levels(x = scedata)

table(scedata@meta.data[["celltype"]])#查看各组细胞数
#BM1  BM2  BM3  GM1  GM2  GM3 
#2754  747 2158 1754 1528 1983
prop.table(table(Idents(scedata)))
table(Idents(scedata), scedata@meta.data[["celltype"]])#各组不同细胞群细胞数
#BM1  BM2  BM3  GM1  GM2  GM3
#  Endothelial  752  244  619  716  906 1084
#  Fibroblast   571  135  520  651  312  286
#  Immune      1220  145  539  270  149  365
#  Epithelial    69   62  286   62   82  113
#  Other        142  161  194   55   79  135
Cellratio <- prop.table(table(Idents(scedata), scedata@meta.data[["celltype"]]), margin = 2)#计算各组样本不同细胞群比例
Cellratio
#BM1        BM2        BM3        GM1        GM2        GM3
#  Endothelial 0.27305737 0.32663989 0.28683967 0.40820981 0.59293194 0.54664650
#  Fibroblast  0.20733479 0.18072289 0.24096386 0.37115165 0.20418848 0.14422592
#  Immune      0.44299201 0.19410977 0.24976830 0.15393387 0.09751309 0.18406455
#  Epithelial  0.02505447 0.08299866 0.13253012 0.03534778 0.05366492 0.05698437
#  Other       0.05156137 0.21552878 0.08989805 0.03135690 0.05170157 0.06807867
Cellratio <- as.data.frame(Cellratio)
colourCount = length(unique(Cellratio$Var1))
library(ggplot2)
pdf("./data/output/细胞比例条形图.pdf",width = 10,height = 6)
ggplot(Cellratio) + 
  geom_bar(aes(x =Var2, y= Freq, fill = Var1),stat = "identity",width = 0.7,size = 0.5,colour = '#222222')+ 
  theme_classic() +
  labs(x='Sample',y = 'Ratio')+
  coord_flip()+
  theme(panel.border = element_rect(fill=NA,color="black", size=0.5, linetype="solid"))
dev.off()


#批量统计图####
table(scedata@meta.data[["celltype"]])#查看各组细胞数
prop.table(table(Idents(scedata)))
table(Idents(scedata), scedata@meta.data[["celltype"]])#各组不同细胞群细胞数
Cellratio <- prop.table(table(Idents(scedata), scedata@meta.data[["celltype"]]), margin = 2)#计算各组样本不同细胞群比例
Cellratio <- data.frame(Cellratio)
library(reshape2)
cellper <- dcast(Cellratio,Var2~Var1, value.var = "Freq")#长数据转为宽数据
rownames(cellper) <- cellper[,1]
cellper <- cellper[,-1]

###添加分组信息
sample <- c("ADJ_ADJ1","ADJ_ADJ2","ADJ_ADJ3","ADJ_ADJ4","ADJ_ADJ5","ADJ_ADJ6",
            "PDAC_PDAC1","PDAC_PDAC2","PDAC_PDAC3","PDAC_PDAC4","PDAC_PDAC5","PDAC_PDAC6")
group <- c("ADJ","ADJ","ADJ","ADJ","ADJ","ADJ",
           "PDAC","PDAC","PDAC","PDAC","PDAC","PDAC")#type
samples <- data.frame(sample, group)#创建数据框

rownames(samples)=samples$sample
cellper$sample <- samples[rownames(cellper),'sample']#R添加列
cellper$group <- samples[rownames(cellper),'group']#R添加列

#调整顺序，把
ADJ <- head(cellper, 6)
PDAC <- slice(cellper, 7:12)
# cellper1 <- bind_rows(ADJ, PDAC)
# cellper2 <- bind_rows(PDAC, ADJ)
cellper <- bind_rows(PDAC, ADJ)

###作图展示
pplist = list()
sce_groups = c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_CTL", "CD8_Trm", "CD8_Te")# 修改
library(ggplot2)
library(dplyr)
library(ggpubr)
library(cowplot)
for(group_ in sce_groups){
  # group_="CD8_CTL"
  
  cellper_  = cellper %>% select(one_of(c('sample','group',group_)))#选择一组数据
  colnames(cellper_) = c('sample','group','percent')#对选择数据列命名
  cellper_$percent = as.numeric(cellper_$percent)#数值型数据
  cellper_ <- cellper_ %>% group_by(group) %>% mutate(upper =  quantile(percent, 0.75), 
                                                      lower = quantile(percent, 0.25),
                                                      mean = mean(percent),
                                                      median = median(percent))#上下分位数
  print(group_)
  print(cellper_$median)
  
  # pp1 = ggplot(cellper_,aes(x=group,y=percent)) + #ggplot作图
  #   geom_jitter(shape = 21,aes(fill=group),width = 0.25) +
  #   stat_summary(fun=mean, geom="point", color="grey60") +
  #   theme_cowplot() +
  #   theme(axis.text = element_text(size = 10),axis.title = element_text(size = 10),legend.text = element_text(size = 10),
  #         legend.title = element_text(size = 10),plot.title = element_text(size = 10,face = 'plain'),legend.position = 'none') +
  #   labs(title = group_,y='Percentage') +
  #   geom_errorbar(aes(ymin = lower, ymax = upper),col = "grey60",width =  1)
  # pp1
  
  pp2 <-
    ggbarplot(
      cellper_,
      x = "group",
      y = "percent",
      add = "mean_se",
      color = "group",
      palette = c("#FC4E07","#00AFBB"),
      # position = position_dodge(1.5),
      xlab = FALSE,
      ylab = FALSE
    ) +
    theme(legend.position = "none") +#不要图例
    labs(title = group_)+#标题
    geom_jitter(shape = 21, aes(fill = group), width = 0.25)
  # +
  #   stat_compare_means(aes(group = "group"), label = "p.signif", label.y = 29)
  # pp2
  
  ###组间t检验分析
  labely = max(cellper_$percent)# 用的max值做的T检验
  compare_means(percent ~ group,  data = cellper_)
  my_comparisons <- list( c("PDAC", "ADJ") )# 修改
  pp3 = pp2 + 
    stat_compare_means(label = "p.format",label.x=2.12)+#size = 3,label.y = 0.6,label.x = 1.7
    stat_compare_means(comparisons = my_comparisons,
                       # label = "p.format",
                       aes(label = paste0("p = ", after_stat(p.format))),
                      # size = 3,
                       method = "wilcox.test",
                       symnum.args = list(cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, Inf), symbols = c("****", "***", "**", "*", "ns")))
  pp3
 
  
  pplist[[group_]] = pp3
}

pdf("./data/output/T细胞批量统计图.pdf",width = 10,height = 12)
# c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_CTL", "CD8_Trm", "CD8_Te")
plot_grid(pplist[['CD4_Tn']],
          pplist[['CD4_Th']],
          pplist[['CD4_Treg']],
          pplist[['CD4_Tem']],
          pplist[['NKT']],
          pplist[['CD8_CTL']],
          pplist[['CD8_Trm']],
          pplist[['CD8_Te']])
dev.off()


